A Detailed Study of Spitzer-IRAC Emission in Herbig-Haro 
Objects (I): Morphology and Flux Ratios of Shocked Emission 



Michihiro Takami 1 , Jennifer L. Karr 1 , Haegon Koh 2 , How-Huan Chen 3 ' 1 , Hsu-Tai Lee 1 



ABSTRACT 

We present a detailed analysis of Spitzer-IRAC images obtained toward six Herbig-Haro 
objects (HH 54/211/212, L 1157/1448, BHR 71). Our analysis includes: (1) comparisons in 
morphology between the four IRAC bands (3.6, 4.5, 5.8 and 8.0 /im), and H2 1-0 S(l) at 2.12 
/im for three out of six objects; (2) measurements of spectral energy distributions (SEDs) at 
selected positions; and (3) comparisons of these results with calculations of thermal H2 emission 
at LTE (207 lines in four bands) and non-LTE (32-45 lines, depending on particle for collisions). 
We show that the morphologies observed at 3.6 and 4.5 /xm are similar to each other, and to 
H2 1-0 S(l). This is well explained by thermal H2 emission at non-LTE if the dissociation rate 
is significantly larger than 0.002-0.02, allowing thermal collisions to be dominated by atomic 
hydrogen. In contrast, the 5.8 and 8.0 /im emission shows different morphologies from the 
others in some regions. This emission appears to be more enhanced at the wakes in bow shocks, 
or less enhanced in patchy structures in the jet. These tendencies are explained by the fact 
that thermal H2 emission in the 5.8 and 8.0 /im band is enhanced in regions at lower densities 
and temperatures. Throughout, the observed similarities and differences in morphology between 
four bands and 1-0 S(l) are well explained by thermal H 2 emission. The observed SEDs are 
categorized into:- (A) those in which the flux monotonically increases with wavelength; and (B) 
those with excess emission at 4.5-/im. The type-A SEDs are explained by thermal H2 emission, in 
particular with simple shock models with a power-law cooling function (A oc T s ). Our calculations 
suggest that the type-B SEDs require extra contaminating emission in the 4.5-/im band. The 
CO vibrational emission is the most promising candidate, and the other contaminants discussed 
to date (H I, [Fe II], fluorescent H2, PAH) are not likely to explain the observed SEDs. 

Subject headings: ISM: Herbig-Haro objects — ISM: line and bands — infrared radiation — ISM: indi- 
vidual (HH 54/211/212, L 1157/1448, BHR 71) 



1. Introduction 

Infrared Array Camera (IRAC) on the Spitzer 
Space Telescope has been extensively used for 
studying young stellar objects (YSOs). The cam- 
era has shown that active YSOs and star form- 
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ing regions are associated with extended infrared 
emission, which is in many cases attributed to 
either shocks (e.g., Noriega-Crespo et al.||2004 



|Smith et a^|2"006l |Walawender et al.||2006j 



Lada 



Teix- 



eira et al.||2008l |Neufeld et alJ|2T)09] |Ybarra and 



120091), scattered continuum in the outflow 
cavity (e.g . |Tobin et al.|2007[|2008| ), or PAH emis- 
sion (e.g., Qiu et al. 2008). Many authors have 



made such identification using three-color images, 
assigning blue, green, and red to 3.6, 4.5, and 8.0 
/im, respectively (see the above references). This 
method is not quantitative, since in principle, indi- 
vidual authors adjust the brightness of each color 
arbitrarily. Even so, stars often appear "blue" in 
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such images due to the fact that the flux is larger 
at shorter wavelengths, while PAH emission ap- 
pears "red" due to excessive emission at 8.0 /xm. 
As a result, the remaining extended emission often 
appears "green". 

The nature of extended infrared emission in 
Spitzer IRAC images has not been fully inves- 
tigated, in particular for this extended "green" 
emission. Shocks often appear green in three-color 
images, and it is often explained by H2 emission 
due to the morphological similarity to H2 1-0 S(l) 
at 2.12 /xm (e .g., |Walawender et al.||2006| |T> 



bin et al. 2007), and some analysis with flux ra- 
tios ( jNeufeld et al.|2009[ |Ybarra and Lada|2009 |. 
In contrast, a number of authors suggest possible 
contribution in shocked emission from additional 
lines and bands, including CO, H I, [Ar II], [Fe 
II], fluorescent H2, and PAH emissions (Noriega- 



Crespo et al.|2004 


Reach et al.||2006| 


Smith et al. 


2006 


Neufeld et al.|2009 


Ybarra and Lada|2009). 



Furthermore, scattered continuum in the outflow 
cavity also tends to appear "green" in three-color 



images (e.g., Tobin et al. 20081. Such an uncer- 
tainty has led to the emergence of a new cate- 
gory of YSOs, namely "extended green objects" 



Cyganowski et al. 2008), whose nature is being 



investigated ( Cyganowski et al.||2009[ ). Better un- 
derstanding of this "green" emission would allow 
further understanding for shock conditions and 
activities for a number of YSOs observed using 
IRAC. 

To investigate the nature of shocked emission 



observed using IRAC, Neufeld et al. (2009) ob- 



served spectra in outflows at 5-37 /xm and dis- 
cussed emission in the four IRAC bands based 
on their calculations of simplified shock models. 
Ybarra and Lada (2009 hearaftcr YL09) have 



made preliminary calculations of molecular hy- 
drogen in non-LTE with isothermal slabs with a 
larger number of transitions, and perform com- 
parisons with observations of HH 54. Follow- 
ing these studies, we present a detailed investiga- 
tion of Spitzer-IRAC emission in shocks associated 
with six Herbig-Haro objects (HH 54/211/212, L 
1157/1448, BHR 71). Our research goal is to, 
with the data set available at present, fully in- 
vestigate the nature of this shocked emission: i.e., 
whether they are explained by thermal H2 emis- 
sion; the implications for temperature (or temper- 
ature structures) and density, and the degree of 



thermalization (LTE or non-LTE) of these regions; 
and possible contaminating emission. The above 
outflows contain a variety of structures such as 
bow shocks, knotty structures, and a bubble-like 
structure, and are thereby ideal for investigating 
emission in shocks for a variety of conditions. 

We investigate morphological similarities and 
differences between the four IRAC bands (3.6, 4.5, 
5.8 and 8.0 /jm), and H2 1-0 S(l) for three out 
of six objects obtained using data obtained by 
Subaru-MOIRCS and VLT-ISAAC. We also clas- 
sify spectral energy distributions (SEDs) observed 
at selected positions. These results are compared 
with:- (1) LTE calculations of thermal H 2 with 
~ 200 lines; and (2) non-LTE calculations of ther- 
mal H2 emission with 32-45 lines for collisions of 
H+He and H2+He. The remaining part of the 
paper is organized as follows. In §2 we describe 
data we have used, including new observations 
using Subaru-MOIRCS. In §3 we show the ob- 
served morphologies, and flux ratios in selected 
positions. In §4 we present our LTE and non- 
LTE calculations, and perform comparisons with 
the observed morphologies and flux ratios. In §5, 
we discuss the nature of shocked emission observed 
in these Herbig-Haro objects, and the applicability 
and limitations of three-color images. 

2. Data 

Archival data in four IRAC bands (3.6, 4.5, 
5.8, and 8.0 /im) were obtained for six Herbig- 
Haro objects: HH 211, HH 212, L 1448, L 1157, 
BHR 71 and HH 54. All the data had been re- 
duced with the post-BCD pipelines developed by 
IPAC. The mean FWHM of the point response 
functions (PRFs) are 1.66, 1.72, 1.88 and 1.98" 
for four bands, respectively. Previous publications 



with these data include Tobin et al. (20071 for L 



1448; Looney et al. (2007) for L1157; Walawen- 



der et al.| ( |2006[ ) for HH 211; |Neufeld et al.| ( |2009[ ) 
for L1448, L1157 and BHR 71; and YL09 for HH 
54. The median of the total integration time per 
pixel ranges 21-322 s among the objects. We did 
not stack all the post-BCD images available on 
archive, but their signal-to-noise are sufficient for 
our study, as shown in later sections. 

Narrow-band images of H2 2.12 /im were ob- 
tained using HH 211, HH 212 and L 1448 us- 



ing MOIRCS (Ichikawa et al. 2006 Suzuki et al 
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2008) on Subaru 8.2-m and ISAAC (iMoorwood 



et al. 19981 on VLT 8.0-m. The MOIRCS obser- 



vations of HH 211 and L 1448 were made on Jan- 
uary 3, 2007. The H 2 filter on MOIRCS provides 
a FHWM wavelength of 0.021 /im. The camera 
consists of two detectors with a FOV of 4'x3' each 
with a pixel sampling of 0".117 pix. -1 , and since 
one of them is sufficient to cover our target, the 
other FOV was used for simultaneously obtaining 
the sky frames. Each object was observed with 
a Detector 1-2-2-1 sequence with an individual 
exposure of 180 s, providing a total on-source in- 
tegration of 720 s. The standard data reduction 
(sky subtraction, flat fielding, field cosmetics, im- 
age stacking) was made using IRAF. Flux calibra- 
tion was made using stars in the FOV and 2MASS 
catalog, with an accuracy of about 10 %. The 1-a 
flux per pixel is 4.1 and 3.6xl0 -9 W m~ 2 str -1 
for HH 211 and L 1448, respectively. The seeing 
of 0".5-0".6 is measured using stars in the same 
FOV. 

The ISAAC images of HH 212 were obtained 



cavity. A similar nebulosity close to the driving 
source is also visible on both sides of L 1157 and 
HH 211. iTobin et al.l ((20071) concluded that the 



through ESO archive (Takami et al. 2006b). These 



images are for two 2'. 5 x 2'. 5 FOV, covering an 
outflow length of about 4'. The exposure time 
of each frame and a median of total on-source 
integration time per pixel are 60 and 840 s, re- 
spectively. The FWHM wavelength of the filter is 
0.028 /im. Flux calibration was made using stars 
in the FOV and 2MASS catalog, with an accuracy 
of about 10 %. The 1-a flux per pixel is 9x 10~ 9 W 
m~ 2 str . The seeing of 0".4 is measured using 
stars in the same FOV. 

3. Results 

3.1. Morphology 

Figure [I] shows the image of the entire outflow 
in the four IRAC bands for L 1448, HH 212, L 1157 
and HH 211. As shown in previous studies, these 
outflows show a collimated morphology with bow 
shocks, knotty, and bubble-like structures. The 
overall morphology in the four bands is similar, 
in particular 3.6 and 4.5 /im. Even so, the mor- 
phology at short (3.6/4.5 /im) and long (5.8/8.0 
/im) wavelengths are not identical as shown later 
in detail. L 1448 exhibits a conical nebulosity just 
north of the driving source, in particular at 3.6 and 
4.5 /im. As discussed by Tobin et al. (20071, it is 



nebulosity at the north end of the L 1448 outflow 
is also due to scattered continuum, by a different 
YSO located there. 

Figure [2] shows the images of BHR 71 and 
HH 54. The BHR 71 outflow shows a collimated 
morphology like the other outflows in Figure [T] 
while the HH 54 outflow exhibits a bright, large 
(~1') bow shock in the northwest and a faint col- 
limated flow toward the southwest. The head of 
the HH 54 bow shock consists of several knotty 
structures, also shown in the 4.5-/iin image in 
YL09. In contrast to the other objects shown in 
Figure [l] the overall morphology of BHR 71 and 
HH 54 is different at short (3.6-4.5 /im) and long 
wavelengths (5.8-8.0 /im). In BHR 71, the short 
wavelengths show a conical nebulosity close to the 
driving source more prominently, while the longer 
wavelengths show the outflow emission in more ex- 
tended regions. Again, the former is presumably 
due to scattered continuum in the outflow cav- 
ity. In HH 54, the emission at shorter wavelengths 
is enhanced at the head of the bow shock, while 
emission at the wakes is also seen at longer wave- 
lengths. 

To show the different morphologies in the four 
bands more clearly, we select several regions and 
show close-up views in Figures [3] and |4j Figure [3] 
shows morphologies of bow shocks at four IRAC 
bands observed in HH 212, BHR 71, and HH 54. 
In HH 212 A, the emission at shorter wavelengths 
is enhanced at the head of the bow shocks, and 
also in a bright wake in the west part of the bow 
shock. In contrast, emission at longer wavelengths 
is enhanced at both wakes of the bow shock. The 
emission at the head of the bow shock is absent at 
8.0 /im. In HH 212 B, BHR 71 A, and HH 54 A, 
the images at all four bands exhibit bright peaks 
at the head of the bow shocks, and the 5.8 and 
8.0 /im images also show faint extended emission 
toward the wakes. In BHR 71 B, the bow shocked 
emission is clearly seen only at 5.8 and 8.0 /im, 
and it is marginal at the other wavelengths. In 
all the regions described above, the wakes of the 
bow shocks are most prominent at 8.0 /im. These 
tendencies are shown in shocked emission with hy- 



drodynamical models by Smith and Rosen ( 2005 1 



explained as scattered continuum in the outflow 



The morphology at 3.6 and 4.5 /im is similar in 
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all of these structures. 

Figure[4]shows the knotty structures seen in HH 
212, L 1448, and L 1157. Emission at 3.6 and 4.5 
/im exhibits patchy structure in these shocks more 
clearly, while it is more blurred at longer wave- 
lengths (5.8-8.0 /im), exhibiting diffuser emission 
components between emission peaks. The ob- 
served morphologies at 3.6 and 4.5 /im are almost 
identical. HH 212 C shows a bow-shock morphol- 
ogy, but the different contrast between the head 
and wakes is not clearly seen, perhaps due to the 
fact that emission is associated with only some 
parts of the wakes. 

Figure[5]shows images of 1-0 S(l), together with 
the IRAC 4.5 and 5.8 /im emission at the selected 
regions in HH 212 and L 1448. To compare the 1-0 
S(l) morphology with the others, we convolve the 
1-0 S(l) images with a gaussian to match the angu- 
lar resolution to IRAC 4.5 /im. These show that 
the morphology at 1-0 S(l) is similar to the 4.5 
/im emission rather than 5.8 /mi. Together with 
Figures [3] and [IJ we conclude that the mophologies 
in 1-0 S(l), IRAC 3.6 /im and 4.5 /im are similar 
to each other. 

3.2. Flux Ratios 

To perform comparisons with model calcula- 
tions in detail, we measured the fluxes at the posi- 
tions shown in Figure[6] These measurements were 
made carefully with the following two caveats. 
First, shocked emission in some objects is severely 
contaminated from diffuse emission in particular 
at 8.0 /zm, presumably due to PAH (see, e.g., 
Reach et al. 20061). Some shocked emission seen 



in at 3.6 and 4.5 /im is contaminated by emission 
from stars, and by scattered continuum near the 
target protostar. We avoided such regions to allow 
reliable discussion in later sections. For each point 
shown in Figure |6j we also carefully subtracted the 
diffuse background emission by measuring it at x- 
and y-directions and fitting them with a 1-D plane. 

Secondly, the different PRFs between the four 
bands would cause systematic errors on flux ra- 
tios. The IRAC Data Handbook issued by Spitzer 
Science Center suggests to perform convolution of 
images to cancel out this effect. We thus convolve 
the 3.6 and 4.5 /im images with PRFs at 5.8 and 
8.0 /im; 5.8 /im images with PRFs at 4.5 and 8.0 
/im; and 8.0 /im images with PRFs at 4.5 and 



5.8 /im. This process yields an effective angular 
resolution of ~4" for each image. Note that the 
PRFs at 3.6 and 4.5 /im are almost identical, thus 
we do not correct the difference between these two 
PRFs. 

Table [I] and Figures [7] show the measured flux 
density at individual positions. Figure [7] shows 
that the observed spectral energy distributions 
(SEDs) fall into two categories: (A) those in 
which the magnitude monotonically decrease as 
the wavelength increases; and (B) those with ex- 
cess emission at 4.5 /im. We define the 4.5 /im ex- 
cess as Am 4 5 = (?Ji 3 6 + m 5 8 ) x 0.5 — m 4 5 , and use 
different grayscales in Figures [7] for different mag- 
nitude of excess. The figure shows that the type- 
B SEDs tend to be brighter than type-A SEDs in 
each. The color gradient of type-A SEDs ranges 
between TO3.6— TO8.o=2.5-5.5, and this corresponds 
the 3.6/8.0-/im flux ratio of 0.03-0.4 based on 
units of MJy str . Thus, the flux monotonically 
increases with wavelength in type-A SEDs. Those 
in L 1157 and HH 54 show that the SEDs in shal- 
low color gradients tend to be brighter than the 
steeper SEDs. 

Figure [8] shows correlations between the ob- 
served IRAC colors and the 4.5-/im excess. Any 
color shows a good correlation with the 4.5-/tm ex- 
cess, implying that the entire SED is characterized 
by this excess. This issue is discussed in detail in 
§5.3. 

Figure [9] show color-color diagrams of these 
SEDs. As expected from Figure [HJ good correla- 
tions are observed between [3.6]-[4.5], [4.5]-[5.8], 
and [5.8]-[8.0] colors. The figure also shows the 
colors of extinction measured in molecular clouds 
by |Chapman et al. ( 2009 ) . Although different ex- 
tinction would produce a color gradient similar 
to type-A SEDs, the slope of the color of type- 
A SEDs in the color-color diagrams in Figure [9] 
is shallower than that expected from extinction. 
Furthermore, the correlation between [4.5]-[5.8] 
and [5.8]-[8.0] can be fitted with a single line with 
both type-A and type-B SEDs. These suggest that 
the observe color gradient is due to different exci- 
tation conditions in shocks, not extinction. 

In addition to the IRAC fluxes, we also mea- 
sured the 1-0 S(l) flux at positions where the data 
are available. Before the flux measuremets, these 
are convolved with IRAC PRFs of 4.5, 5.8 and 8.0 
/im to approximately match the angular resolu- 
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tions. The results are shown in Table [2] Figure 10 
shows the IRAC SEDs at these positions normal- 
ized to 1-0 S(l) emission. The 4.5 /xm/1-0 S(l) flux 
ratio is higher for SEDs with a larger 4.5-^tm ex- 
cess ATO4.5. More quantitatively, the flux ratio in 
SEDs with AWI4.5 > 0.3 is larger than those with 
A7714.5 < 0.3, typically by a factor of ~2. 

4. Comparisons with Modeled H2 Emis- 
sion 

We compare the results in §3 with calcula- 
tions of thermal H2 emission. These calcula- 
tions include:- (1) LTE calculations with 245 en- 
ergy levels, including 207 lines in the four IRAC 
bands (§4.1); and (2) non-LTE calculations with 
36-49 levels, including 32-45 lines in the four 
IRAC bands (§4.2). The non-LTE calculations 
are made for collisional particles for two extreme 
cases: H+He and H 2 +He. For the temperature 
structure, we use simplified shock models with a 
power-law cooling rate (A oc T s ), in addition to 
isothermal cases. 

The comparison between observations are made 
with SEDs and color-color diagrams. Although 
the flux ratio with 1-0 S(l) could be useful for 



testing models (Figure 10), we do not perform 



such comparisons due to uncertainties in extinc- 
tion, and the absolute flux calibration of the IRAC 
data for extended sources (see IRAC Data Hand- 
book). 

4.1. LTE models 

4- 1.1. Description and basic characteristics 

The list of energy levels and transitions are ob- 
tained from models of photo-dissociation regions 



by Draine and Bcrtoldi ( 1996 ). This covers 245 en- 



ergy levels at .E/&:=0-43000 K, and among them, 
we use all of the transitions at 3.09-10.2 /im (207 
lines) . For the ^-coefficients we use recent calcula- 
tions by Wolniewicz et al. ( 1998 ) . For isothermal 



cases, the IRAC flux is calculated from each spec- 
trum with the following equation: 



I IRAC — £l ^ Rl(kjk) 

j,k 



A 



■3 k- 



n H2 ASl gje 
Z{T) 



(1) 



where Iirac is the flux observed in a IRAC band; 
Cj is the conversion factor from electrons to MJy 
str -1 (sec IRAC Data Handbook); j and k are the 
upper and lower levels of transitions, respectively; 
-Rj(Ajfc) is the spectral response function at the 
wavelength Ajk', Ajk is the Einstein A-coefficient; 
nn 2 is the column density of molecular hydrogen; 
A is the area of the telescope mirror (0.4636 m 2 ); 

is the angular area corresponding to a pixel 
(1". 22x1". 22 or 3.498X10" 11 str); gj is the sta- 
tistical weight of the level j; Ej is the upper level 
energy; and Z(T) is the partition function. 

In addition to the isothermal cases, we also 
calculate the IRAC fluxes assuming the temper- 



ature structures used in|Brand et al. (1988| here- 
after B88), and |Neufeld and Yuan| ([2008[ here- 



after NY08). B88 developed a simplified dissocia- 
tive J-shock model with a power-law cooling rate 
A = AqT s to explain infrared H 2 emission at a 
broad range of wavelengths. The column density 
at the energy level j is given as: 

Nj/gj cx Je- T ^ T [Z(T)Aer)]- 1 dT, 

cx J e~ T i /T Z{TY X T~ s dT, (2) 

where Tj = Ej/k. The authors provide an approx- 
imate formula for the population of each energy 
level of:- 



N s /g s 



t; 



(T 1 +r v y 



(3) 



with T„=6000 K. The B88 model has been suc- 
cessful in roughly explaining near-infrared spectra 



of H 2 emission in shocks (e.g., Brand et al. 



1988 



Everett et ah] [19951 |Richter et al.||1995[ |Takami| 
et al.||2006a| ) The power index s is often close to 
4.7, i.e., cooling dominated by molecular hydrogen 
(Brand et al. 1988). This model can also explain 



the relatively constant temperature measured at 
a variety of shock regions with near-infrared lines 



(T ~ 2000 K, e.g., Burton et al. 1989 


Eisloffel 


et al.|2000 


Nisini et al. 2002; Takami et a 


1. 2006a; 


Beck et al. 


2008 


). The NY08 model is essentially 



the same as B88, but Equation (2) is numeri- 
cally integrated with the upper limit of tempera- 
ture T=4000 K. This temperature approximately 
correspond to the dissociation temperature of H 2 
molecules (e.g., Lepp and McCray 1983). This 



model has been successfully applied to H 2 spectra 
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at the 5-28 (im range in some shocks, including 
the supernova remnant IC 443 (NY08) and several 
Herbig-Haro objects QNeufeld et al.|2009[ ). We set 
the lowest temperature of 30 K for numerical in- 
tegration. 

Figure [□] shows the SEDs of LTE H 2 for 
isothermal cases, and B88 and NY08 models with 
different power indexes for the cooling function. 
For B88 and NY08, the flux monotonically changes 
with wavelength, exhibiting a tendency similar to 
the type-A SEDs described in §3.2. The slope of 
the SEDs is steeper with a larger power index s 
because of the contribution of gas at low temper- 
atures to longer wavelengths. The B88 and NY08 
models provide almost the same results at large s, 
while the SED is steeper for NY08 with small s, 
due to a smaller contribution of high-temperature 
slabs to the latter. None of the LTE models dis- 
cussed here explains the type-B SEDs, which show 
excess emission at 4.5 /im. 



Figure 12 shows the IRAC and 1-0 S(l) fluxes as 
a function of temperature. These are normalized 
at 2000 K, i.e., the typical temperature observed 
in shocks based on ground-based infrared spec- 



troscopy ( 


e.g., Burton et al. 1989 Eisloffel et al. 


2000 


Nisini et al.|2002| Takami et al.|2006a[ Beck 


et al. 2008 


). The figure shows that, compared with 



1-0 S(l), the 3.6 and 4.5 /xm emission enhances 
regions at higher temperatures; the 8-/xm emis- 
sion does regions at lower temperatures; and the 
5.8-/im emission traces regions similar to 1-0 S(l). 
This contrasts with the observed tendency for the 
3.6 and 4.5 /im emission to show morphologies 
similar to 1-0 S(l), rather than 5.8-/im (§3.1). 

4-1-2. Color-color diagrams 



Figure 13 shows color-color diagrams of ob- 
served and modeled fluxes. In the [4.5]-[5.8] vs. 
[3.6]-[4.5] diagrams, the models could explain the 
observed colors for only a fraction of SEDs ([4.5]- 
[5.8]=0.4-l.l). Isothermal models at 2000-3000 
K and the NY08 models with s ~ 3.0 explain the 
observed flux at [4.5]-[5.8]~0.5. The B88 models 
explain the observed flux ratio at [4. 5] - [5. 8] ~ 1.0 if 
we assume a large extinction (Ak=5). In any case, 
the model curves cannot entirely fit the distribu- 
tion of observed colors in the diagrams. In par- 
ticular, the observed fluxes show [3.6]— [4.5] > 1.3, 
while any LTE models flux show [3.6]-[4.5]<1.3 if 
extinction is not included. 



Regarding the [5.8]-[8.0] vs. [4.5]-[5.8] dia- 
grams, the B88 and NY08 models explain the color 
of the observed SEDs with little or no 4.5 /im ex- 
cess (Am4. 5 < 0.3) well. The SEDs at a larger 4.5 
/im excess show smaller [4.5]-[5.8] colors by up to 
~1.0 mag. The isothermal models cannot explain 
either type-A or type-B SEDs. 

4.2. Non-LTE models 

4-2.1. Description and Characteristics 

Following YL09, we adopt A-coeffcients pro- 
vided by Wolniewicz et al. (19981, and collisional 



rate coefficients with H 2 and He by |Le Bourlot 



et al. (1999). For collisional rate coefficients with 



H, we adopt Wrathmall et al. (2007) and Le Bour- 



lot et al. (1999), and show the results separately. 
According to Wrathmall et al. (2007), they pro- 



vide rate coefficients with a better accuracy than 



Le Bourlot et al. ( 1999 ) because of improved rep- 
resentation of the vibration eigen functions. In 
contrast, coefficients provided by |Le Bourlot et ah] 
(1999) can explain the observed SEDs better as 



shown in §4.2.2. Table [3] summarizes the energy 
levels included in the online tables of IWrathmalll 
et al.l d2007b andlLe Bourlot et a.1.1 (11999b . Consid- 



ering the completeness of the transitions covered 
for ortho- and para-H 2 , we include the lowest 49 
levels (E/k = 2 x 10 4 ) in case of H+He collisions, 
and 36 levels (E/k = 1.69 x 10 4 K) in case of 
H2+He collisions. For the Hc/H abundance ratio, 
we adopt 0.098 based on Dappen (2000). 



Transitions and rate coefficients between ortho- 
and para-H 2 are not included in the above phys- 
ical parameters, hence we calculate the popula- 
tions of ortho- and para-H 2 separately, and com- 
bine them assuming an ortho/para ratio of 3. This 
ratio is corroborated by near-infrared spectroscopy 
of shocked regions (see, e.g 



Everett et al. 



Giannini et Tab] [2001] |Nisini et al.||2002| |Takarm| 



et al. 2006a Beck et al. 2008), and explained by 



1995 



chemical processes of grain surfaces sufficient for 
thermal equilibrium ( Le Bourlot et al.||1999 ). 

The calculations are made for two cases: for 
isothermal and the NY08 temperature structure 
(§4.1.1). Due to the limited number of transitions 
included, non-LTE calculations produce large er- 
rors at large temperatures (or small s) at 3.6 /im. 
To investigate this, we compare our non-LTE cal- 
culations at the highest density regime (i.e., at 
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LTE) with that of §4.1, i.e., calculations with all 



the transitions used in Draine and Bertoldi ( 1996 1 



We then define limitations of parameter ranges 
to guarantee a calculation error within ^20 %. 
These are: an upper limit of T of 2000-3000 K 
for isothermal models; and a lower limit of s of 3.0 
and 4.7 for NY08 models with 49 and 36 energy 
levels, respectively (Appendix A). The errors for 
the other wavelengths (4.5, 5.8, and 8.0 /im) are 
much smaller, allowing us to model the emission 
with sufficient parameter ranges for our interest. 



Figure 14 shows the IRAC and 1-0 S(l) fluxes 
as a function of hydrogen number density for three 
cases: (1) collisions with H+He, adopting the col- 

( [1999] ; 



lisional rate coefficients of Lc Bourlot et al. 



(2) same as (1) but adopting the collisional rate 



coefficient for H of Wrathmall et al. (20071; and 



(3) collisions with H2+He , adopting the collisional 
rate coefficients of Le Bourlot et al. (1999). These 
calculations are made at temperature T=2000 K. 
For all cases the emission at shorter wavelengths 
requires higher densities for thermalization. The 
5.8 /im emission reaches LTE at higher density 
than 8.0 /un; 3.6 and 4.5 /im emission reaches 
LTE at even higher density. The 3.6 /jm emission 
reaches thermal equilibrium at a much higher den- 
sity than 4.5 /xm emission for H2+He collisions, 
while the difference between two lines are rather 
marginal for H+He collisions. The dependence of 
the 3.6 (im emission on density is almost identical 
with 1-0 S(l). The non-LTE models with H+Hc 
collisions best explain the the observed morpho- 
logical similarities between 3.6 /im, 4.5 /jm, and 
1-0 S(l) emission shown in §3.1. 

The above discussions are often made with criti- 
cal densities, although there are several difficulties 
for applying this to IRAC emission. These are:- 
(1) each band includes multiple lines with differ- 
ent critical densities; and (2) ro-vibrationally ex- 
cited states have two critical densities, representa- 
tive of the thermalization of rotational and vibra- 
tional states, respectively. We thus alternatively 
define the critical density of IRAC band emission 
as I{n cr it) — 0.5 x Ilte- This definition is based 
on the fact that the equation yields the correct 
critical density for the 2-level model. 

Table [4] shows the critical densities of H2 emis- 
sion for four IRAC bands and 1-0 S(l) based on 
this definition (T=2000 K). As expected from Fig- 
ure [141 the densities are significantly higher at 



shorter wavelengths: the 5.8 /im emission having 
a larger critical density than 8.0 /im by a factor 
of 6-11; the 4.5 /im having a higher critical den- 
sity than the 5.8 /im by a factor of 4-11. The 
critical densities for 3.6 and 4.5 /jm are similar 
to each other in case of H+He collisions, while the 
former is larger than the latter by a factor of 14- 
15 in case of H 2 +He collisions. In addition to the 
H+He and H 2 +He gas, we also show the critical 
densities without collisions with He. These show 
that collisions with He are not dominant compared 
with those with He or H2, but not negligible, in 
particular for H 2 +He gas (20-30 %). Two differ- 



ent tables for H collisions ( Le Bourlot et al7||1999 
Wrathmall et al.|[2007 ) yield critical densities dif- 
ferent by a factor of 2-10. 

The critical densities for H+He and H 2 +He col- 
lisions are significantly different, particularly at 
shorter wavelengths. The former is larger than the 
latter by by a factor of 50-500 for 3.6 /im, and 3-6 
for 8.0 /im. This implies that the 3.6 /im flux is 
determined by collisions with H even with a low H 2 
dissociation rate (>0.002-0.02). The 8.0 /im flux 
is determined by collisions with H in cases the dis- 
sociation rate is much larger than 0.2-0.4. 



Figure[15]shows the SEDs of non-LTE and LTE 
gas with different particles and collisional rate 
coefficients. Since the shorter wavelength IRAC 
bands require higher densities for thermalization, 
these are fainter in the SEDs at low densities. As 
clearly seen in the figure, the deficit at lower densi- 
ties is largest at 4.5 /xm for collisions with H+He, 
while it is largest at 3.6 /im for collisions with 
H 2 +He. In any case, the SEDs do not show 4.5 
/im excess observed in type-B SEDs. 

4.2.2. Color-color diagrams 

Figure[l6]shows color-color diagrams for the ob- 
served SEDs, and modeled emission for isothermal 
cases. The [4.5]-[5.8] vs [3.6]-[4.5] diagrams show 
that the colors of SEDs with weak/no 4.5-/im ex- 
cess are entirely explained either with (1) H+He 
collisions based on collisional rate coefficients by 
Le Bourlot et al.| ( |1999) , or (2) H 2 +Hc collisions. 
In the former case, the different colors observed 
in these SEDs (Am.4.5 < 0.3) are attributed to 
different densities but with a single temperature 
(T ~ 2000 K). In the latter case, for H 2 +He col- 
lisions, the variation of the observed colors are 
attributed to a combination of different density 
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and temperature in the case of H2+He collisions. 
These models indicate a lower limit of hydrogen 
number density of ^10 5 cm -3 for H+He collisions, 
and ~10 7 cm~ 3 for H2+H2 collisions. The H col- 



lisions based on parameters by |Wrathmall et al. 



(2007) explain only a fraction of the SEDs ob- 
served ([4.5]-[5.8]=0.0-1.2). None of the non-LTE 
models with an isothermal temperature structure 
explain the observed colors plotted in the [5.8]- 
[8.0] vs [4.5]-[5.8] diagrams. 

Figure [IT] shows color-color diagrams for the 
observed SEDs, and modeled emission based on 
NY08. As in the isothermal cases, the [4.5]-[5.8] vs 
[3.6]-[4.5] diagrams show that the colors of SEDs 
with weak/no 4.5-/im excess are explained with 
either H+He collisions based on |Le Bourlot et al.| 
(1999), or H2+He collisions. In both cases, the 



variation of the observed colors is attributed to a 
combination of different densities and power in- 
dexes s. In the [5.8]-[8.0] vs. [4.5]-[5.8] diagrams, 
the observed colors are well fitted by the models 
at LTE. This is in contrast to the fact that the 
observations are explained by non-LTE regimes in 
the [4.5]-[5.8] vs [3.6]-[4.5] diagrams. This dis- 
crepancy is attributed to different critical densi- 
ties between the four bands. As for the isothermal 
models, the [4.5]-[5.8] vs [3.6]-[4.5] diagrams indi- 
cates a lower limit for the hydrogen number den- 
sity of ~10 5 cm~ 3 for H+Hc collisions, and ~10 7 
cm -3 for H 2 +He collisions. Such densities yield 
close to LTE values in the [5.8]-[8.0] vs. [4.5]-[5.8] 
diagrams. 

5. Discussion 

We summarize and discuss how the observed 
morphologies and SEDs are explained with models 
of H2 emission, and implications for the physical 
conditions. Discussions of the morphologies and 
two types of SEDs are shown in §5.1-5.3, respec- 
tively. In §5.4, we discuss the applicability and 
limitations of thee-color images, which have been 
used by many authors. 

5.1. Observed Morphologies vs. H2 Mod- 
els 

As shown in §3.1, the morphologies observed at 
3.6 /im, 4.5 /im and 1-0 S(l) match well, while 
the 5.8 and 8.0 /im emission show slightly differ- 
ent morphologies. In regions where we can clearly 



see the morphological differences, the former is as- 
sociated with the head of a bow shock or they 
exhibit clumpy/pachy structures; the latter is en- 
hanced at the wakes of bow shocks or they exhibit 
more blurred structures. As discussed in §4.1.1 
and §4.2.1, the similar morphologies of 3.6 /im, 
4.5 /im and 1-0 S(l) emission is best explained by 
non-LTE models with H+He collisions. Since the 
collisional rate coefficients for H are much larger 
than H2 , it would require only a small dissociation 
rate. As discussed in §4.2.1, a dissociation rate 
of 3>0. 002-0. 02 would allow collisional excitation 
and de-excitation dominated by H at T=2000 K, 
i.e., a typical temperature observed in Figure [T6| 

The morphological differences between 4.5, 5.8, 
and 8.0 /im emission are attributed to either dif- 
ferent temperatures or densities. In bow shocks, 
higher excitation is expected at the head than the 
wakes due to different shock velocities (e.g. 



Harti- 



gan et al.|[T987| |Morse et al.1|i"993l |Allen and B~ 



ton 1993). The morphological differences in IRAC 



emission observed in such shocks are attributed to 
the fact that the thermal H2 emission at shorter 
wavelengths tends to be associated with regions at 
higher temperatures (§4.1.1, §4.2.2). The fact that 
the emission at shorter wavelengths is enhanced in 
patchy structures may be due to different densi- 
ties, rather than temperatures. This agrees with 
our calculations, which show that emissions at 
shorter wavelengths have higher critical densities 
(§4.2.1). 

5.2. Type-A SEDs vs. H 2 Models 

As shown in §3.2, the observed SEDs are cat- 
egorized into: (A) those in which the flux mono- 
tonically increases with wavelength; and (B) those 
with excess emission at 4.5 /im. As discussed 
in §4.2.2, SEDs with no or weak 4.5-/xm excess 
(type-A or weak type-B) are explained by H 2 emis- 
sion, in particular with simple shock models with 
a power-law cooling (A oc T s ). Either collisions 
with H+He or H2+He can explain the observed 
colors, but the former required the collisional rate 



coefficients of H of Le Bourlot et al. ( 1999 ) , rather 
than more recent calculations of Wrathmall et al.l 
(2007). The color-color diagrams show that the 



IRAC emissions at short wavelengths (in particu- 
lar 3.6 /im) originate from H2 in non-LTE, while 
that at long wavelengths (in particular 5.8 and 8.0 
/im) is due to H2 at LTE. This discrepancy is at- 
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tributed to different critical densities between four 
IRAC bands (§4.2.2). 

Wrathmall et al.| ( |2007 ) state that they provide 



improved results compared with |Le Bourlot et al. 
(1999). However, models with their coefficients 



explain only a fraction of observed IRAC colors 
at short wavelengths (§4.2.2 and Figure 17). If 
their coefficients are accurate, the observed SEDs 
may be explained as follows. Thermal collisions 
are dominated by H2+HC for SEDs with steep 
gradients ([4.5]-[5.8]>1.2), but are dominated by 
H(+He) for SEDs with shallow gradients ([4.5]- 
[5.8]<1.2). This explanation is consistent with the 
fact that a higher dissociation rate is expected in 
regions at higher temperatures. 



If we assume H2+He collisions, Figure 17 



plies a gas density significantly higher than that 



provided by Neufeld et al. (20091, who estimated 



the density in several outflows with low- J lines as- 
suming H2+He collisions. In Figure [FT] observed 
SEDs with weak or no 4.5-/zm excess (ATO4.5 < 
0.3) are explained with a density of more than 10 7 
cm" 3 . In contrast, Neufeld et al. (2009) fit their 
spectra with much smaller densities (~10 4 cm" 3 ). 
This discrepancy is consistent with the fact that 
our selection of positions is biased to relatively 
bright regions in all four IRAC bands. Due to a 
high critical density at short wavelengths (10 6 -10 7 
cm -3 for 3.6 and 4.5 /im emission for H 2 -|-He col- 
lisions), such positions should be associated with 
gas at high densities. This selection bias would 
also explain why a different power-law index s is 
required to explain our SEDs (3.0-6.0) compared 
to the low- J spectra observed by |Neufeld et al.| 
( [2009] ) (2-3). 

The temperature structure we used (A cx T s 
at T max =4000 K) corresponds to a simplified J- 
shock, although our calculations would not reject 
C-shocks for explaining the type-A SEDs. Indeed, 
the isothermal models also account for the ob- 
served colors in the [4.5]-[5.8] vs. [3.6]-[4.5] di- 
agrams, as shown in Figures 13 and [17] This 
suggests that the temperature structure near the 
highest temperature, which affects the colors at 
short wavelengths, is arbitrary to some extent. 
Thus, a flat-top temperature structure predicted 



by C-shock models (e.g., 


Draine et al.|1983 


Kauf- 


man and Neufeld 1996 


Le Bourlot et al. 


2002) 



may also fit well the type-A SEDs. Testing this 
scenario is beyond the scope of this paper. 



5.3. Type-B SEDs 

Our results suggest that SEDs with a large 4.5- 
/im excess (strong type-B) cannot be solely at- 
tributed to thermal H2 emission. Modeled SEDs 
at either LTE and non-LTE do not show such 
excess emission (§4.1.1 and §4.2.1). The [5.8]- 
[8.0] vs. [4.5]-[5.8] diagrams in Figures 13 and 



[17] show discrepancies between the observed SEDs 
with ATO4.5 > 0.3 and models up to [4.5]-[5.8]~l. 
These discrepancies between the strong type-B 
SEDs and H 2 models are in contrast to the analysis 
by YL09, who explain such colors in shocks associ- 
ated with HH 54 ([3.6]-[4.5] >1.0, [4.5]-[5.8] <0.7) 
with non-LTE H2 calculations with the lowest 49 
levels. However, the temperature YL09 measured 
is 2000-4000 K, implying that the [3.6]-[4.5] color 
they model would significantly suffer from errors 
due to the limited number of transitions (up to 
~0.6 mag, see Appendix A). 

Thus, we suggest a different origin, or the 
presence of an additional emission component for 
the type-B SEDs. The continuous color varia- 
tion shown in Figure [8] suggests continuous varia- 
tion in density and temperature in the observed 
regions. The contamination of CO vibrational 



emission, suggested by several authors (Noriega 



Crespo et aL][2004l |Reach et alT]|2006| |Qiu et al 
2008), is the most likely candidate to explain 
all the observed tendencies. The entire type- 
B SEDs, which show (1) excess emission at 4.5- 
/im; and (2) shallower magnitude gradient at the 
other wavelengths (Figures [7] and [Sj, are explained 
with a combination of density and temperature 
higher than type-A SEDs. The former is explained 
with the fact that the CO emission contributes to 



mostly the 4.5-/im flux (see Reach et al. 20061, 
and that the contribution is comparable to H2 
emission if the number density of H atoms is 10 7 - 
10 s cm" 3 (see Figure 11 of NY08). The latter is 
explained by the fact that the SEDs for thermal H2 
emission show shallower gradients in magnitude at 
a high density or temperature. 

Alternative candidates for the 4.5-fim excess 
emission could be H I (in particular Br a) and [Fe 



II] (|Noriega-Crespo et al.|2004| [Reach et al.||2006 



Qiu et al. 2008 YL09) . Both emission mechanisms 



are expected in dissociative shocks (e.g., Hollen- 
bach and McKee|[l989| . |Hollenbach and McKee 
(1989) show that such shocks diminish thermal 
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H2 emission at T > 1000 K, since reformation of 
most of the H2 molecules occurs at lower temper- 
atures. This trend is clearly seen in bow shocks in 



the Orion-KL outflow by Tedds et al. (19991, who 



show that the [Fe II] 1.64 /im emission is associ- 
ated with the head, while the H2 2.12 /im emis- 
sion is associated with the wakes. Such a deficit 
of the observing flux is not observed in the type-B 
SEDs: these are comparable to or brighter in all of 
the four IRAC bands than type- A SEDs. We thus 
conclude that the excess emission at 4.5-/im is not 
due to either H I or [Fe II]. 

Based on IRAC observations, Smith et al.| 



(2006) reported the presence of UV fluorescent 
H2 emission in shocks associated with the DR 
21 outflow. To investigate this possibility, we 
calculated the flux of photodissocoiation regions 



(PDRs) modeled by Draine and Bertoldi (1996) 
with a variety of densities (n# = 10^ to 10 b cm -3 ) 
and UV field (\ — 1 to 10 5 ). These arc shown in 
the color-color diagrams in Figure [Tsj together 
with the observed colors in shocks. The observed 
colors show no correlation with the PDR models, 
implying that contribution of UV fluorescent H 2 
emission to our SEDs is negligible. 

5.4. Discriminating Shocks from PAH and 
Stars with IRAC Colors 

As described in §1, three-color images (blue, 
green, red for 3.6, 4.5 and 8.0 /im, respectively) 
have been used by a number of authors to discrim- 
inate between stars, PAH, shocks and scattered 
continuum in star forming regions. Stars are rel- 
atively bright at shorter wavelengths hence these 
appear blue in such images; PAH emission has a 
remarkable excess at 8.0 /im hence it appears red; 
and the others tend to appear green. 

Figure [19] shows the 3.6, 4.5, and 8.0 /im col- 
ors observed in shocks, PAH in NGC 7023 esti- 
mated using ISO spectra (IRAC Data Handbook), 
and stars. The figure shows that some shocks can 
look "red" in the three-color images. As shown 
in §4, such shocks have relatively low tempera- 
ture and/or density. The figure also shows that 
the [3.6]-[4.5] color allows us to clearly discrimi- 
nate shocks ([3.6]-[4.5]>l.l) from stars and PAH 
([3.6]-[4.5]<l.l). The observed [3.6]-[4.5] color is 
lowest for the SEDs with [4.5]-[8.0] =1.0-2.0, and 
larger [3.6]-[4.5] colors for the other SEDs are due 
to either their low densities (§4.2.2) or contami- 



nating emission at 4.5 /im (§5.1). The presence 
of a lower limit in the [3.6]-[4.5] color is also sup- 
ported by empirical use of 4.5/3. 6-/im flux ratio to 
identify shocks in some recent studies (Qiu et al. 
20081 |Zhang and Wang||2009 ). 
Figure 
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shows the three-color image and 
4.5/3.6-/im flux ratio in the southern end of the 
L 1448 jet. The brightness scale of the three-color 
image is automatically adjusted using ds9 devel- 
oped by Smithsonian Astrophysical Observatory. 
The figure shows that this shock region contain 
both greenish and reddish components, agreeing 
with the above discussion with Figure [19] The 
flux ratio map shows that the ratio is comparable 
or higher than ~ 1.5 (corresponding to [3.6]-[4.5] 
> 0.9) in the entire jet, in contrast to the stars in 
the same image (flux ratio <C 1.5). 

YL09 state the possibility of PAH emission in 
shocks seen in IRAC images. Our results do not 
show clear evidence for such contaminating emis- 
sion. In Figure [l9j the colors measured in shocks 
are distributed perpendicularly to the color of 
PAH, with a significant color offset. 

6. Summary and Conclusions 

We present detailed analysis of shocked emis- 
sion in four Spitzer-IRAC bands associated with 
six Herbig-Haro objects (HH 212, L1448, L1157, 
HH 211, BHR 71 and HH 54). Our research goal 
is to fully investigate the nature of this shocked 
emission: i.e., is this explained by emission from 
molecular hydrogen; what are the implications for 
temperature and density of these regions; what 
are the possible sources of contaminating emission. 
We compare the observed morphologies and flux 
ratios measured at high-S/N regions with models 
of thermally excited molecular hydrogen emission 
with different excitation conditions and two types 
of temperature structures: isothermal and power- 
law cooling (A cx T s ). For HH 212, L 1448 and HH 
211, narrow-band images of 1-0 S(l) observed us- 
ing Subaru-MOIRCS and VLT-ISAAC are added 
for our study. 

The morphologies observed at 3.6 /im and 4.5 
/im and 1-0 S(l) match well, while the 5.8 and 
8.0 /im emission show different morphologies. In 
regions where we can clearly see the morphological 
differences, the former is either associated with the 
head of bow shock or they exhibit clumpy /pachy 
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structures; the latter is either enhanced at the 
wakes of bow shocks or they exhibit more blurred 
structures. The observed SEDs in selected regions 
are categorized into: (A) those in which the mag- 
nitude monotonically increases with wavelength; 
and (B) similar to (A) but with excess emission 
at 4.5 /mi. The type-B SEDs tend to be brighter 
than type-A in each object, and the magnitude of 
excess is well correlated with the color of the entire 
SED. 

We compare the above results with (1) LTE 
calculations with 245 energy levels, including 207 
lines in four IRAC bands; and (2) non-LTE calcu- 
lations with 36-49 levels, including 32-45 lines in 
four IRAC bands. The non-LTE calculations are 
made for collisional particles for two cases (H+Hc 
and H2+He). These show that H2 emission ex- 
plains the morphological similarities between 3.6 
fj,ra and 4.5 /im and 1-0 S(l), and their differ- 
ences from 5.8 and 8.0 /im emission. The for- 
mer fact is particularly well explained with non- 
LTE H2 if the dissociation rate is much larger than 
0.002-0.02, allowing thermal collisions dominated 
by atomic hydrogen. The latter fact is explained 
by the fact that H2 emission at longer wavelength 
bands should be enhanced at lower temperatures 
and/or densities. 

The H2 models also successfully explain type- 
A SEDs. The models with power-law cooling 
rate explain different colors between SEDs, with 
different densities and power indexes particularly 
well. These are explained by collisions with either 
H2+He, or H+He with collisional rate coefficients 
of H provided by ILe Bourlot et all ((19991. The 



variation in SEDs shows that the colors at short 
wavelengths are attributed to non-LTE H2, while 
those at long wavelengths are attributed to LTE 
H2 due to their lower critical densities. In contrast 
to the type-A SEDs, our calculations suggest that 
thermal H 2 cannot explain the 4.5-/xm excess in 
type-B SEDs. We suggest that the combination of 
CO and H2 emission qualitatively explains the ob- 
served tendencies. It is not likely that H I, [Fe II], 
fluorescent H2 or PAH significantly contributes to 
the SEDs we observed. 

Three-color images (blue, green, red for 3.6, 4.5 
and 8.0 /im, respectively) are widely used by many 
authors to instantly discriminate between shocks, 
PAH emission and stars. We show that shocked 
emission at low temperatures (i.e., large power- 



law index for cooling function) can be "red" in 
such images, analogous to PAH emission. Our 
new and thorough measurements of colors suggest 
that these can be discriminated with the [3. 6] -[4. 5] 
color: >1.0 for shocks and <1.0 for PAH. 
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A. Error of non-LTE calculations with limited number of transitions 

Our non-LTE calculations of thermal H 2 include the lowest 49 levels (E u /k < 2 x 10 4 K) for thermal 
collisions with H or H+He , and the lowest 36 levels (E u /k < 1.69 x 10 4 K) for collisions with H2 or H2+He. 
In practice, the transitions from higher energy levels contribute to the total IRAC flux, thus our non-LTE 
calculations would have larger errors in particular at short wavelengths and high temperatures. 

To investigate the missing flux from the upper levels, we calculate the fraction of the LTE flux from 
the lowest 36-49 levels, and compare them with the LTE flux we derive in §4.1.1. Figure [21] shows the 
difference of SEDs in these cases, at different temperatures and power indexes for cooling. Tables [5] shows 
the contribution of the lowest 36-49 levels to the total flux at different temperatures, and different power 
indexes for the cooling function (s). These show that the 3.6-/zm flux would significantly suffer from errors 
due to limited numbers of transitions. The lowest 36-49 levels are esponsible for only 20-50 % of the total 
flux at T=4000 K; 57-80 % for s=3.0. Table [6] shows the estimated errors for IRAC colors in magnitudes. 
The parameter range of non-LTE calculations in Figures [TB] and [17| are determined based on these tables. 

Note that the population of the lower 36-49 levels is still responsible for a large fraction of population of 
H2 (94-97 % at T=4000 K). The significant contribution from the upper-level transitions to the 3.6-/zm flux 
is due to a large number of lines (~100) and their large A coefficients. 
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Fig. 1 — IRAC Images of L1448, HH 212, L 1157, and HH 211. The fl ux is shown for each image with units 
of MJy str -1 . The cross indicates the position of the driving source ( [Lee et al.||2008| |2009| |Bachiller et al. 
2001 Girart and Acord||2001 ). Boxes with dashed lines indicate the regions shown in Figures |3j|5j Archival 
data is missing in the south part (0'.5-l') of the L1448 jet at 4.5 /im. 
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Fig. 2. — Same as Fig ure [T] but for BHR 71 and HH 54. The position of the driving sources of BHR 71 are 
based on ( Chen et al. 2008 ) . The driving source of the HH 54 outflow is not clear (see Ybarra and Lada 
20091). 
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Fig. 3. Bow shocks and similar shock structures observed in HH 212, BHR 71 and HH 54. All the figures 
are rotated to align the direction of the down stream to the left. 
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Fig. 4. — Clumpy shock structures at selected regions in HH 212, L 1448 and L 1157 B. All the figures are 
rotated to align the direction of the down stream to the left. Archival data is missing near the right end of 
the L1448 B at 4.5 fim. 
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Fig. 6. — Positions where we measure spectral energy distributions in IRAC bands. See text for details. 
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Fig. 7. — Flux observed in Figure [6j Different grayscale and dots are used for different degree of excess 
emission at 4.5 /im. The excess ATO4.5 is defined as (m3.6+m5.8)x0.5-TO4.5. The error bars are shown only 
for those larger than the size of the dots. 
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Fig. 8. — Correlation between the observed IR color and excess at 4.5 /xm. The error bars are shown only 
for those larger than the size of the dots. 
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Fig. 9. — Color-color diagrams of the observed flux. The black and gray dots indicate the results with excess 
4.5-^m magnitude Ato 4 . 5 > 0.3 and < 0.3, respectively. The arrows for extinction is based on Chapman et 
al. (2009) for A K > 2. 
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Fig. 10. — Observed IRAC flux normalized by the 1-0 S(l) flux. The left and right figures show SEDs with 
different excess 4.5-^m magnitudes (Am4. 5 > 0.4 and < 0.4, respectively). Different grayscale and dots are 
used for different degree of excess emission at 4.5 /im. 
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Fig. 11. — Spectral energy distributions for LTE models per unit column density of H2. (left) For isothermal 
cases at T=1000-4000 K. (middle) For the Brand ct al. (1988) models with different power indexes of the 
cooling function. Normalized by the column density of excited H2, since the models yields an infinite column 
density for the ground state (see Equation 3). (right) Same as the middle but for the Neufeld & Yuan (2008) 
models. 
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Fig. 12.— IRAC and 1-0 S(l) fluxes for LTE-H 2 
as a function of temperature. The flux for each 
band or line is normalized to that at T=2000 K. 
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Fig. 13. — Same as Figure [9] but the colors of LTE models are overplotted with a solid curve/line and square 
dots, (left) At a single temperature, (middle) For the Brand et al. (1988) models with different power 
indexes of the cooling function. The gray dashed line with square dots indicate the case for A#-=5. (right) 
For Neufeld & Yuan (2008) models with different power indexes of the cooling function. 
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Fig. 14. — IRAC flux of non-LTE H2 at T=2000 K as a function of hydrogen number density, (left) For 
collisions with H+He, with the collisional rate coefficients of Lc Boutlot et al. (1999) for both particles; 
(middle) Same as the left figure but the collisional rate coefficients of Wrathmall et al. (2007) are used for 
collisions with H; (right) For collisions with H 2 and He. All the fluxes are normalized to those at LTE. 
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Fig. 15.— Spectral energy distributions of non-LTE and LTE H 2 in IRAC bands at T=2000 K. (left) For 
collisions with H+He, with the collisional rate coefficients if Le Boutlot et al. (1999) for both particles; 
(middle) Same as the left figure but the collisional rate coefficients of Wrathmall et al. (2007) are used for 
collisions with H; (right) For collisions with H 2 and He. SEDs for non-LTE cases are shown for hydrogen 
number density of 10 3 , 10 4 , and 10 5 cm~ 3 . 
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Fig. 16. — Color-color diagrams for observations and non-LTE models for isothermal cases, (left) For 
collisions with H+He, with the collisional rate coefficients by Le Boutlot et al. (1999) for both particles; 
(middle) Same as the left figure but the collisional rate coefficients by Wrathmall et al. (2007) are used for 
collisions with H; (right) For collisions with H2 and He. The upper and lower plots are diagram with [4.5]- 
[5.8] vs. [3.6]-[4.5] and [5.8]-[8.0] vs. [4.5]-[5.8], respectively. In the upper-left and upper-middle figures, 
the modeled colors for T=3000 K are shown in dashed line, since it would have a relatively large error for 
calculations (>0.2 mag. for [3.6]-[4.5], see Appendix A). 
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Fig. 17. — Same as Figure but with temperature structures used in Neufeld & Yuan (2008). Solid and dotted 
lines in the upper figure shows the modeled values with constant power indexes s and densities, respectively. 
In the lower figures, the solid lines are used for models with constant densities to clearly show that the LTE 
values match the observations. In the upper-right figure modeled colors at s—3.8 is shown in dashed lines 
due to a relatively large uncertainty (>0.2 mag. for [3.6]-[4.5], see Appendix A). 
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Fig. 18. — Color-color diagrams for SEDs observed in shocks (black and gray dots), and PDR models 
provided by Draine & Bertoldi (1996) at n H = 10 2 to 10 6 cm~ 3 and x = 1 to 10 5 (squares). 
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Fig. 19. — 3.6-, 4.5-, and 8.0-/xm colors of shocks, stars, and PAH. The corresponding color for the three- 
color images are also shown beside the labels for the horizontal and vertical axes. Extinction added to the 
stars is based on Chapman et al. (2009) for A K > 2. Color of PAH in NGC 7023 is based on IRAC Data 
Handbook. 
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Fig. 20. — (left) Three-color image at the south 
end of L 1448. The color scale is automatically ad- 
justed using ds9, developed by Smithsonian Astro- 
physical Observatory, (right) The 4.5/3. 6-/im flux 
ratio. The measurements are made at positions 
where the signal-to noise exceeds 10-ct both at 3.6 
and 4.5 /mi. Before taking the flux ratio, the 3.6- 
/im image is convolved with the PRF at 4.5 /im, 
and 4.5-/im image is convolved with the PRF at 
3.6 /im, to develop higher signal-to-noise. 
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rig. 21. — IRAC flux of LTE H 2 calculated with the lowest 36/49 levels, and all the transitions included 
in Draine & Bertoldi (1996). (upper left) For the isothermal case at 7=2000 K; (upper right) same but 
7=4000 K; (lower left) Neufeld & Yuan (2008) model with s=4.7; (lower right) same but s=3.0. 
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18. UD 
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16.63 
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4 
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1D.D I 


( ^ft f\l\ 
[<V.Uo ) 


10. Oo 


( ^ft ftQ^ 

^<U.Uo ) 


L 1157 


1 


17.52 


(<0.03) 


16.10 


(<0.03) 


15.68 


(<0.03) 


14.87 


(<0.03) 





z 


18. L I 


^<U.Uo ) 


1 7 ft/i 

1 ( .U4 


(<0.03) 


16.23 


{ ^ft ftQ^ 
^<U.Uo J 


10. oO 


( ^ft 

l^<U.Uo J 




Q 




id orr 


( ft ft/1 ^ 
^ U.U4 J 


10 in 
18. 1U 


(<0.03) 


17.10 


{ ^ft ftQ^ 
^<U.Uo J 


1 R Oft 
1D./U 


^ ^ft ft^ 
l^<U.Uo J 




A 

4 


m 01 

ly.oi 


/ ^ A ft o\ 

(<U.v6 ) 


18. oz 


(<0.03) 


16.98 


/ ft^ 
(<U.Uo J 


10. Dy 


^ft ftQ\ 
(<U.UO J 




c 



on 1 k 


( ri ft.7 ^ 

^ U.U ( ) 


1 Q cn 

18. oy 


(<0.03) 


17.31 


1 ^ft ftQ^ 
(^<U.Uo J 


1 K ftft 

10. yy 


^^ft ftQ^ 

^<U.Uo ) 


HH 211 


1 


16.80 


(<0.03) 


14.62 


(<0.03) 


14.93 


(<0.03) 


14.68 


(<0.03) 




2 


18.08 


( 0.03 ) 


16.92 


( 0.04 ) 


16.41 


( 0.04 ) 


15.53 


( 0.06 ) 




3 


17.53 


( 0.07 ) 


15.80 


(<0.03) 


15.91 


( 0.05 ) 


15.58 


( 0.12 ) 




4 


17.77 


( 0.08 ) 


16.38 


( 0.04 ) 


16.18 


( 0.05 ) 


15.65 


( 0.10 ) 


BHR 71 


1 


18.88 


( 0.15 ) 


17.63 


( 0.08 ) 


16.86 


( 0.05 ) 


16.23 


( 0.07 ) 




2 


17.69 


( 0.05 ) 


16.56 


(<0.03) 


15.74 


(<0.03) 


15.11 


(<0.03) 




3 


16.63 


(<0.03) 


14.95 


(<0.03) 


15.22 


(<0.03) 


15.09 


(<0.03) 




4 


17.68 


(<0.03) 


15.38 


(<0.03) 


15.80 


(<0.03) 


15.46 


(<0.03) 




5 


17.43 


( 0.06 ) 


15.97 


(<0.03) 


15.44 


(<0.03) 


14.71 


(<0.03) 




6 


17.68 


( 0.08 ) 


16.10 


(<0.03) 


15.55 


(<0.03) 


14.81 


(<0.03) 




7 


17.44 


(<0.03) 


15.99 


(<0.03) 


16.00 


(<0.03) 


15.56 


( 0.04 ) 


HH 54 


1 


17.72 


(<0.03) 


16.39 


(<0.03) 


15.30 


(<0.03) 


14.21 


(<0.03) 




2 


17.36 


(<0.03) 


16.10 


(<0.03) 


15.03 


(<0.03) 


14.10 


(<0.03) 




3 


18.01 


( 0.03 ) 


16.44 


(<0.03) 


15.21 


(<0.03) 


14.33 


(<0.03) 




4 


17.84 


(<0.03) 


16.31 


(<0.03) 


15.32 


(<0.03) 


14.66 


(<0.03) 




5 


17.71 


(<0.03) 


16.12 


(<0.03) 


15.62 


(<0.03) 


15.01 


(<0.03) 




6 


20.52 


( 0.35 ) 


18.73 


( 0.08 ) 


16.99 


( 0.04 ) 


15.57 


(<0.03) 




7 


20.34 


( 0.29 ) 


18.64 


( 0.07 ) 


17.03 


( 0.04 ) 


15.65 


(<0.03) 




8 


20.96 


( 0.37 ) 


19.16 


( 0.06 ) 


17.64 


( 0.04 ) 


16.60 


(<0.03) 




9 


18.81 


( 0.09 ) 


17.49 


( 0.03 ) 


16.76 


( 0.03 ) 


16.11 


( 0.04 ) 




10 


19.40 


( 0.16 ) 


18.01 


( 0.05 ) 


16.79 


( 0.03 ) 


16.17 


( 0.04 ) 
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Table 2: Observed 1-0 S(l) Flux 



Object Position Flux 
(W in' 2 str' 1 ) 



HH212 


1 


l.lxlCT 6 




2 


7.1xlCT 7 




3 


2.ixicr 6 




4 


2.7xlCT 7 


L1448 


1 


2.9xlCT 6 




2 


4.1xlCT 6 




3 


l.OxlO" 6 


HH211 


1 


7.3xl(T 6 




2 


1.6xlCT 6 




3 


3.3xlCT 6 




4 


2.7xlCT 6 



Table 3: Summary of Transitions included in Le Bourlot ct al. (1999) and Wrathmall et al. (2007) 



Paper 


Particle 


Counterpart 
for collisions 


Total energy 
levels included 


Completeness of Transitions 


Lo Bourlot (1999) 


ortho-H2 


H, He. para-R~2 


24 


up to 23 levels (E u /k < 19086 K) 






ortho-H2 


22 


up to 22 levels (E u /k < 18979 K) 




para-H~2 


H, He , para-H2 


27 


up to 26 levels (E u /k < 19112 K) 






ortho-H2 


19 


up to 19 levels (E u /k < 16880 K) 


Wrathmall ot al. (2007) 


ortho-R~2 


H 


54 


up to 54 levels (E u /k < 34037 K) 


para-H 2 


H 


54 


up to 54 levels (E u /k < 32127 K) 



Table 4: Critical hydrogen nucleus densityin cm 3 defined by I /Ilte=0-5 (T=2000 K) 



Collisions with 


Rcfcrcncc a 




IRAC Band 




1-0 S(l) 






[3.6] 


[4.5] 


[5.8] 


[8.0] 


H 


W07 


2.3 x 10 4 


3.6 x 10 4 


8.7 x 10 3 


1.5 x 10 3 


1.7 x 10 4 


H+He 


W07 (H) + L99 (He) 


2.3 x 10 4 


3.6 x 10 4 


8.5 x 10 3 


1.4 x 10 3 


1.7 x 10 4 


H 


L99 


2.2 x 10 5 


1.5 x 10 5 


2.3 x 10 4 


3.4 x 10 3 


1.6 x 10 5 


H+Hc 


L99 


2.2 x 10 5 


1.4 x 10 5 


2.1 x 10 4 


3.0 x 10 3 


1.6 x 10 5 


H 2 


L99 


1.8 x 10 7 


1.2 x 10 6 


1.1 x 10 5 


1.0 x 10 4 


1.7 x 10 7 


H 2 +He 


L99 


1.2 x 10 7 


8.8 x 10 5 


8.0 x 10 4 


7.8 x 10 3 


1.1 x 10 7 



a W07 ... Wrathmall et al. (2007); L99 ... Le Bourlot et al. (1999) 



Table 5: Contribution of the lowest 36 and 49 levels to the total LTE flux (per cent) 

Lowest 36 Levels Lowest 49 Levels 

(E u /k < 16900 K) (E u /k < 20000 K) 





[3.6] 


[4.5] 


[5.8] 


[8.0] 


[3.6] 


[4.5] 


[5.8] 


[8.0] 


T(K) 


















1000 


99.9 


99.9 


>99.9 


>99.9 


>99.9 


>99.9 


>99.9 


>99.9 


2000 


84 


97 


99.5 


>99.0 


96 


99.2 


99.8 


>99.9 


3000 


44 


89 


96 


99.6 


75 


95 


98.0 


99.6 


4000 


21 


78 


89 


98.5 


53 


88 


94 


98.5 


power index s 


















3.0 


57 


92 


98.8 


>99.9 


80 


96 


99.3 


>99.9 


3.8 


65 


93 


99.2 


>99.9 


84 


98 


99.5 


>99.9 


4.7 


75 


94 


99.5 


>99.9 


89 


97.5 


99.7 


>99.9 


6.0 


87 


96 


99.8 


>99.9 


95 


98.4 


99.9 


>99.9 
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Table 6: Error of IRAC Color for LTE Calculations with the lowest 36/49 Levels 





Lowest 36 Levels 


Lowest 49 Levels 




(E u /k < 16900 K) 


(E u /k < 20000 K) 




[3.6]-[4.5] [4.5]-[5.8] [5.8]-[8.0] 


[3.6]-[4.5] [4.5]-[5.8] [5.8]-[8.0] 


T(K) 






1000 


<0.01 <0.01 <0.01 


<0.01 <0.01 <0.01 


2000 


0.15 0.03 0.01 


0.03 0.01 <0.01 


3000 


0.76 0.09 0.04 


0.26 0.03 0.02 


4000 


1.43 0.15 0.11 


0.56 0.07 0.06 


power index s 






3.0 


0.51 0.08 0.02 


0.20 0.03 0.01 


3.8 


0.38 0.07 0.01 


0.15 0.03 0.01 


4.7 


0.25 0.06 0.01 


0.10 0.02 <0.01 


6.0 


0.10 0.04 <0.01 


0.04 0.02 <0.01 
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